Cyclization of a Polymer: A First Passage Problem for a Non-Markovian Process 
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We discuss a problem of cyclization of a polymer molecule, which is an important example of 
reaction in a system showing strongly non-Markovian behavior on the timescales of interest. We 
show that the knowledge of the joint three-time probability distribution of the end-to-end distance 
is sufficient for the full description of the cyclization kinetics, so that the survival probability follows 
rigorously as a solution of the Volterra integral equation. The corresponding kinetics can easily be 
evaluated numerically. We moreover discuss how do some well-known approximations appear from 
this exact scheme due to decoupling. 
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Kinetics of reactions involving polymers (especially the 
luminescent energy transfer) has attracted much theo- 
retical attention in the last decade (see e.g. ||]), due 
to experimental relevance and to a plenitude of applica- 
tions. Cyclization of a polymer, or luminescence quench- 
ing in a system of a donor and acceptor attached to the 
ends of the polymer molecule belongs to the most com- 
plex problems in the field: It is a case of a diffusion- 
controlled reaction in a system exhibiting considerably 
non-Markovian behavior on the timescales of interest 
i, |, |, |, |, 0, i- This non-Markovian behavior mir- 
rors the essentially many-particle nature of the system. 
Even for the simplest situation, in which the polymer is 
modelled by a Rouse chain, i.e. the excluded volume ef- 
fects and the hydrodynamical interactions between the 
monomers ("beads") are neglected, the existing theoreti- 
cal approaches fail to give an exact description of the sit- 
uation and make additional physical assumptions. Thus, 
their domain of applicability is inevitably restricted. In 
what follows we give an exact solution for a special sit- 
uation when the reaction takes place on a contact, and 
show that it only depends on three-time joint probability 
density of the end-to-end distance. We also discuss how 
do some standard approximations emerge. 

The theoretical approaches typically start from the 
Langevin equation for the whole chain oi N ^ 1 
monomers, and try to reduce the overall dynamics to the 
manageable effective one. The simplest approach consists 
in describing the evolution of the end-to-end distance r 
as an effective diffusion JH. A more accurate Wilemski- 
Fixman theory (WF) [^[ and its equivalents D ex- 
press the reaction rates through the effective sink-sink 
correlation function and assume local equilibration. The 
comparison of theoretical results with the results of nu- 
merical simulations of the underlying dynamics for 
the case when the reaction takes place with probability 
1 whenever the ends of the chain (the monomers and 
A'') approach each other at the distance a, shows that 
the WF approach leads to systematic overestimate of the 
mean first contact time of the ends of the polymer, while 
neglecting the effective diffusion assumption leads to un- 
derestimates. The simulations in Q started from the 
equilibrium (Boltzmann) initial distribution of the end- 



to-end distance tq, with rg > a. Another approach ||^ is 
based on the reduction of the full equation of motion to 
the equation of motion of the end-to-end distance using 
the projection operator formalism. The corresponding 
equation is again approximate. We note that the meth- 
ods of Refs. [|[ ^, ^, H are based on the Smoluchowski 
approach in the theory of diffusion-controlled reactions, 
in which the reaction probability follows from the solu- 
tion of the full or reduced diffusion equation. 

In the case of the reaction taking place with probabil- 
ity 1 on contact, as discussed in Ref. 0, introducing a 
sink term in the exact or reduced Fokker-Planck equa- 
tion (like in the WF theory or in Ref. ||]) is equivalent 
to posing an absorbing boundary condition P(r — a) = 
(a " black sphere" ) . Such problem can be reformulated as 
the first passage problem from the initial distance r = rg 
between the ends to r = a In the Markovian case, 

such first passage problem can be solved using an alterna- 
tive approach |^ based on the renewal property; the idea 
stemming from MontroU and Weiss This approach 
can be generalized to the non-Markovian situation. 

In what follows we concentrate on the solution of the 
initial-condition problem, not averaged over the initial 
positions. The reason is two-fold: First, the recent suc- 
cess of micromanipulation of polymers makes the situ- 
ation experimentally relevant, while no theoretical pre- 
dictions for the kinetic curves are available. Second, the 
initial condition problem renders clear the geometrical 
nature of standard approximations 1^, |4[ ^ . The results 
for the equilibrium initial condition can be obtained by 
additional numerical averaging. 

Let us first note that the reaction process depends only 
on the single variable r so that the problem is essentially 
one-dimensional. Moreover, the trajectories of our pro- 
cess are continuous and nowhere differentiable, just like 
the trajectories of the usual Wiener process are {vide in- 
fra). Let us consider the relation between the distribu- 
tion of the first passage time to a sphere of the radius 
r = a around the origin of the coordinates F{a, t \ rg, 0) 
and the conditional probability for the ends to be found 
at distance r: If the trajectory starts at point rg, and is 
found at r at time t, it may have already visited r before, 
at some time t' < t. Thus, the conditional probability to 
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be at r at time t, provided the particle started at ro at 
t = 0, G(r, t I To, io)? is given by the foUowing equation: 



G(r,i I rQ,0) = S{r-ro)S{t) 



(1) 



F{r,t' I ro,0)G(r,i I r, t'; ro, 0)dt', 



where G{r,t \ r, t';ro,0) is the conditional probability to 
be at r at time t, provided r was visited earlier at time t' 
and that the particle started at at t = 0. This equa- 
tion is essentially the definition of the first passage time 
distribution and holds for all processes with continuous 
trajectories, whether Markovian or not. For a Marko- 
vian process the conditional probability G{r,t \ ro,0) is 
a Green's function, and G{r,t \ r,t';rQ,0) depends only 
on the latest the arguments to the right of the line, so 
that G{r,t \ r,t';ro,0) = G{r,t | r,t'). In this case, our 
Eq.(|l]) reduces to a well-known renewal equation for the 
first-passage time Q. If the particle definitely docs not 
start at rg = r, the (5-functional term can be omitted. 

According to the Bayes formula, G{r,t \ ro,0) = 
P(r,t;ro,0)/P(ro,0) and G{r,t \ r,t';ro,0) - 
P(r,t;r,i';ro,0)/P(r',t';ro,0) where P(ro,0), 
P{r,t;rQ,0) and P{r,t;r,t';rQ,0) are the one-, two- 
and tree-time joint probability distributions. The 
corresponding one-dimensional distributions as functions 
of r can be expressed as integrals of the overall joint 
probabilities over the surface of the sphere of radius 
r, so that P{r,t \ ro,to) = (isP(r, i; ro, to) and 
P(r, t; r, t'; ro, to) = /g (isds'P(r, t; r, t'; ro, to), so that 
the overall equation reads 



/^dsP(r,t;ro,0) 



^'(ro,0) 



(2) 



P(r,t' I ro,0) 



Jg Jg dsds'P{r, t; r, t'; ro, 0) 
/5ds'P(r',t';ro,0) 



dt'. 



(where |ro| > r). It is reasonable to rewrite this equation 
in the following form: 



/ P(r,t'|ro,0)g(t,t',ro)dt' = l, 
Jo 



(3) 



with the kernel 



Q(t,t',ro) = 



Jg Jg dsds'P(r, t; r, t'; ro, 0)P(ro, 0) 



[Jg ds'P(r', t'; ro, 0)] [Jg rfsP(r, t; ro, 0)] 



(4) 

This is an exact equation expressing the pdf of first pas- 
sage time through the three-time joint probability distri- 
bution of the end-to-end distance. As stated, we concen- 
trate on the initial value problem and discuss the distri- 
bution of the first passage times as depending on ro. 

Note that for the Rouse chain and many other so-called 
generalized Gaussian structures Q the random process 
describing temporal changes of the distance between each 
two beads is Gaussian, since it is a weighted sum of many 



Gaussian random variables describing the uncorrelated 
displacements of beads. Thus, for example, the three- 
time joint probability distribution is 



P(r, t; r, t'; ro, 0) - ±^ exp 
Vdet A 



^RA 



(5) 



where the vector R is a 9-component vector 
{x,y,z,x',y',z',XQ,yQ,zo) and the covariance matrix 
A = [{RiRj)] consists of 9 diagonal blocks: 



D(0) D(t-t') D(t) 
A = I D(t - t') D(0) D(t') 
D(t) D(t') D(0) 



(6) 



where D(t) is a diagonal matrix D(t) — i4>{t) (I is a 
unit matrix). The function (j){t) is a relaxation function 
of the structure, (/)(t) = {x{t)x{0)). The corresponding 
two-point distributions read 



with 



P(r,t;ro,0) 



B 



(2^)- 



= exp 



VdetB 



D(0) D(t) 
D(t) D(0) 



^RB ^R 



(7) 



(8) 



and the one-point density reads P(ro,0) = 

[27rmr'^'exp[-{rl/2m)]- 

The integrals over the two-point functions, say 
Jg dsP{r,t;ro,0), can be evaluated analytically. Taking 
the z-axis to follow the direction of ro one gets (r=a): 



/ dsP(r,t;ro,0) = 
Js 



(9) 



IS 

2„2 exp 



"2[0^(O)-0^(t)] 



sinh 



aro4>(t) 



TT 



aro(/.(t)V02(O)-02(t) 



Let us now consider the three-time pdf 
Jg Jg dsds' P{r,t;r,t' -jTojO). Taking the point ro to 
lay on the z-axis and the end-point r to have a zero 
y-coordinate, the integrals over the azimuthal angles 
can be performed analytically, so that the corresponding 
distribution reduces to 

^4 



P(r, t;r,t';ro,to) = 



(27r)5/2A3/2{^} 



(10) 



exp 



a2(202_02_03)+^2(^2_^2) 



Iq 



2A{^} 

a? sin Tp s\nd{4>o4>i - 4>24'^) 



A{^} 



X exp [/(a, ro, {(?!)})/A(?/))] sinV" sin^JdV'rf^?, 

where d and ■0 are the corresponding polar angles. Here 

A{(j)} = 0|] -I- 201 02 03 -"/"o (01+ 02 +'^'3), -^0 is the modified 
Bessel function, and 

/(a,ro,{0}) = cos 01 cos i?(0203 - 0o0i) + 

-l-aro [(0102 - 0003) cos -h (0103 - 0O02) cos-!/(]ll) 
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where (po = 0(0), 0i = (j){t — t'), (f)2 = 4>{t) and (j)^ = (p{t'). 
These forms are universal and apply to any Gaussian 
model, whether linear or branched chain, or a network. 

Let us now turn to a special case of the Rouse chain. 
For long enough chain the time correlation function 
can be approximated through 



= —J2p'^ exp{-Xpt) 



(12) 



where for a Rouse chain the summation runs over the 
odd integers p and Ap = P^/tr In what fol- 

1 /2 

lows we use the rms size of the chain (i^) and the 
Rouse time as length and time units, so that r and 
t are nondimensional. In these units we get 0(t) = 
■§j E^r=o(2" + 1) exp(-(2n + This function is con- 

tinuous in zero, so that (^[x{t) — x{0)]'^'^ = 2[l — (f>{t)] 

tends to zero for small t, which fact verifies the conti- 
nuity of the trajectories necessary for our consideration. 
The Volterra equation can be then solved numerically 
by approximating the integral in Eq.(^ by a finite sum. 
This solution and the evaluation of the double integral, 
Eq.(|lO|) in the kernel is performed using MATHCAD. The 
corresponding first-passage time distributions for differ- 
ent initial conditions are shown in Fig.l for the initial 
distances rg = 1,2,4 and 8 and for a = 0.1. The numer- 
ical accuracy of the results is around one per cent. We 
note that with decreasing the distance, the typical mean 
first passage time decreases very fast; the maximum of 
the curve corresponding to rg — 1 is at <,nax ~ 0.01, and 
is not resolved on the scales of Fig.l. This effect is due 
to the small-scale motion corresponding to higher modes 
iQ, ^ and mirrors the compact exploration of space by 
the chain's end. On the other hand, the large-scale re- 
laxation is mostly connected with the lowest mode, which 
is slow. 

Let us now discuss the emergence of the WF- 
approximation and the Markovian approximation of Ref. 
1^ within our scheme and return to Eq.(||). The WF 
approximation corresponds to a pseudo-Markovian de- 
coupling G{r,t I r,i';ro,0) ~ G{r,t \ r,t') in Eq.(0), 
which, in the language of P is: P(r, t; r, i'; ro, 0) ~ 
P(r, t- r, t')P[vo,0) and P(r', i'; Tq, 0) ^ P(r', t')P{vo, 0). 
This corresponds to the approximation of the kernel Q 
by 



}wF{t, t',rQ) = 



JgJgdsds'Pir,t;r,t')P{ro,0) 
[J^ds'P{r',t')] [JgdsP{r,t;ro,0) 



(13) 

Note that this approximation has a purely geometrical 
nature: It assumes that the distance from the point Tq to 
all points on the surface S is approximately the same and 
corresponds to using the mean value theorem in evaluat- 
ing the integrals in the r.h.s. of Eq.(|^). This assumption 
is exact for tq = and for tq — > oo; in our case rg > a 
it delivers a low-order expression in a/rQ. The solutions 




FIG. 1: The first passage time pdf's for a = 0.1 and vq = 
1, 2, 4 and 8 (from left to right). The solid curves correspond 
to the solution of equation, Eq.(^ with the exact kernel, 
Eq.(Q). The dotted lines (shown for ro = 1. 2. 4') correspond 
to the Wilemski-Fixman approximation, Eq.(|l3|), see text for 
details. The inset shows the case a = 0.5 and ro = 2, for 
which the WF-approximation ceases to perform reasonably. 



of Eq.(g) with the kernel Qwf are shown in Fig.l for 
ro = 1, 2 and 4 with dashed lines. We note that for a/ro 
small the approximation performs excellently and leads 
to slight overestimate of the first passage times. The 
discrepancies grow when a/ro gets larger: fixing for ex- 
ample, r = 2 and increasing the value of to a = 0.5 we 
find that the height of the maximum is underestimated 
by 15% on the cost of the somewhat fatter tail of the 
distribution; the mean value of the first passage time is 
overestimated by around 3%. For even larger a/ro the 
accuracy of the approximation worsens rapidly. 

We now show how do the usual equations of the 
Wilemski-Fixmann theory appear from this decoupling 
approximation. To do this we return to Eq.(^, now read- 
ing as 



JgdsP{r,t;ro,0) 



F{r,t' I ro,0) 



Jg Jgdsds'P{r,t; r, t'] 
Jgds'P{r',t') 



(14) 



dt'. 



and average both sides over the equilibrium distribution 
of ro, P(ro,0). We note that the function multiplying 
F under the integral in the r.h.s. does not depend on 
ro, so that only F is averaged, giving rise to the av- 
eraged first-passage distribution F, and that the aver- 
age in the l.h.s. is a one-time marginal distribution, i.e. 
Jg dsP{r, t). Thus, the equation for F reads: 



Jg Jgdsds'P{r,t]r,t') 



^^"'''^'''°^JsdsPir,t)Jgds^P{r',t^) 



dt' = 1. 



(15) 

Noting that the process r(t) is stationary we denote 
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Jg Jg dsds' P{r,t;r,t') = C{t — t'). Due to correlation 
decoupling at long time, Jg dsP{r, t) ds'P{r' , t') is ex- 
actly C(oo) — Coo, if the corresponding limit does not 
vanish. Thus, we can write 



' F{T,t'\x„Q)^^^-^dt' 
C [oo) 



1. 



(16) 



Applying the Laplace transform to the both sides of the 
equation, we get 1/u = F{u)C{u)/Coo, so that 



F{u) = Coo/ [uC{u)] 



(17) 



Moreover, since Coo — linit^oo C'(t) , for u small one 
has C{u) — Coo/u + A + with A — hm„^o(C('u) — 
Coo/u). Noting that the mean first passage time r = 



J^-'tF{t)dt^-^F{u) 



one arrives at 



u=0 




A 



(18) 



which is exactly Eq. (12) of Ref. . The approach of Ref. 

is an approximation of the same nature: It assumes 
that the random process r(t) is an Ornstein-Uhlenbeck 
process (the only one Gaussian Markovian process, the 
one with (f>{t) ~ exp(— at), with a = Q ), for which 
Eq.(|l4|) is exact. Compared to the WF-theory it contains 
an additional assumption, namely one of the effectively 
exponential relaxation of and has typically a lower ac- 
curacy. 

Compared to the general case, our solution applies 
to a rather special sink function (" black sphere" , or a 
d -sink), but our derivation here is rigorous and much 
less technical than the original WF approximation, and 
can apply to a variety of cases where the approximation 
fails. Let us say a few words on the general sink func- 
tion. In general, in the WF theory one has C{t — t') — 
Jy Jydrdr'S{r)S{r')P{r,t;r,t'), where S{r) character- 
izes the sink strength. Although such expressions have no 
immediate counterparts within the first passage time for- 
malism, they can be obtained as mean-field approxima- 
tions when one associates S{r) with the density of absorb- 
ing regions in a problem where the absorbing boundary 
is corrugated or even not singly-connected. In this case, 
the corresponding surface integrals Eq.(jl^) arc changed 
for the volume integrals weighted with the sink function. 



The same approximation can be done, of course, also 
in a general non-Markovian form, Eqs.(^jJ). The corre- 
sponding forms may be no more exact, but still take into 
account mode details of the process than WF ones do. 

Let us summarize our findings. The cyclization of a 
polymer is an example of a diffusion-controlled reaction 
in a system exhibiting strongly non-Markovian behav- 
ior on the timescales of interest. Considering the reac- 
tion as a first-passage process, one is able to provide an 
exact equation governing the reaction probability. This 
equation can be readily solved numerically, giving the ki- 
netic curves. We also discuss how the popular Wilemski- 
Fixman approximation emerges from the exact scheme 
as a pseudo-Markovian decoupling. 
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